# Copyright (c) HySoP 2011-2024
#
# This file is part of HySoP software.
# See "https://particle_methods.gricad-pages.univ-grenoble-alpes.fr/hysop-doc/"
# for further info.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
"""
@file stretching_dir.py
Directional stretching frontend (operator generator).
"""
import sympy as sm
from hysop.constants import DirectionLabels, Implementation, StretchingFormulation
from hysop.tools.htypes import check_instance, to_tuple, to_list, first_not_None
from hysop.tools.decorators import debug
from hysop.tools.numpywrappers import npw
from hysop.parameters.scalar_parameter import ScalarParameter
from hysop.core.graph.graph import generated
from hysop.fields.continuous_field import Field
from hysop.topology.topology import Topology
from hysop.topology.cartesian_descriptor import CartesianTopologyDescriptors
from hysop.operator.directional.symbolic_dir import DirectionalSymbolic
from hysop.symbolic.relational import Assignment
from hysop.operator.directional.directional import DirectionalOperatorFrontend
[docs]
class DirectionalStretching(DirectionalSymbolic):
"""
Directional stretching using the symbolic code generation framework.
"""
@debug
def __new__(
cls,
formulation,
velocity,
vorticity,
variables,
dt,
C=None,
A=None,
name=None,
implementation=None,
base_kwds=None,
**kwds,
):
base_kwds = first_not_None(base_kwds, {})
return super().__new__(
cls,
name=name,
variables=variables,
dt=dt,
candidate_input_tensors=None,
candidate_output_tensors=None,
base_kwds=base_kwds,
implementation=implementation,
exprs=None,
**kwds,
)
@debug
def __init__(
self,
formulation,
velocity,
vorticity,
variables,
dt,
C=None,
A=None,
name=None,
implementation=None,
base_kwds=None,
**kwds,
):
"""
Initialize directional stretching operator frontend.
Vortex stretching is the lengthening of vortices in three-dimensional fluid flow,
associated with a corresponding increase of the component of vorticity in the stretching
direction due to the conservation of angular momentum.
Solves dW/dt = C * f(U,W)
where U is the velocity, W the vorticity and f is one of the following formulation:
CONSERVATIVE FORMULATION: f(U,W) = div( U : W )
MIXED GRADIENT FORMULATION: f(U,W) = (A*grad(U) + (1-A)*grad(U)^T) . W
where : represents the outer product U:W = (Ui*Wj)ij.
. represents the matrix-vector product.
^T is the transposition operator
U and W are always three dimensional fields.
C is scalar or 3d vector-like of symbolic coefficients.
A is a scalar symbolic coefficient.
f(U,W) is always directionally splittable.
Parameters
----------
formulation: hysop.constants.StretchingFormulation
The formulation of this stretching operator:
CONSERVATIVE: f(U,W) = div( U : W )
GRAD_UW: f(U,W) = grad(U) . W
GRAD_UW_T: f(U,W) = grad(U)^T . W
MIXED_GRAD_UW: f(U,W) = (A*grad(U) + (1-A)*grad(U)^T) . W
where A is a user given scalar or 3d vector like that contains
scalar (floating point) coefficient(s) or parameter(s).
The most common mixed gradient formulation is A=0.5.
GRAD_UW, GRAD_UW_T and MIXED_GRAD_UW formulations are only valid
for incompressible flows, see Notes.
velocity: Field
Velocity U as read-only input three-dimensional continuous field.
vorticity: Field
Vorticity W as read-write three-dimensional continuous field.
variables: dict
Dictionary of fields as keys and topology descriptors as values.
dt: ScalarParameter
Timestep parameter that will be used for time integration.
C: scalar or vector like of 3 symbolic coefficients, optional
The stretching leading coefficient C can be scalar or vector like.
Contained values should be numerical coefficients, parameters or
generic symbolic expressions such that C*f(U,W) is directionally splittable.
Here * is the classical elementwise multiplication.
Default value is 1.
A: scalar symbolic coefficient, optional
Should only be given for MIXED_GRAD_UW formulations.
ValueError will be raised on other formulations.
The linear combination coefficients A is a scalar.
Contained value should be a numerical coefficient, a parameter (or a generic symbolic
expression) such that C*(A*grad(U).W + (1-A)*grad^T(U).W) is directionally splittable.
Here * is the classical elementwise multiplication.
Default value is 0.5.
name: str, optional, defaults to 'stretching'.
Name of this stretching operator.
implementation: Implementation, optional, defaults to None
Target implementation, should be contained in available_implementations().
If None, implementation will be set to default_implementation().
base_kwds: dict, optional, defaults to None
Base class keywords arguments.
If None, an empty dict will be passed.
kwds:
Keywords arguments that will be passed towards implementation
operator __init__.
Notes
-----
The MIXED GRADIENT formulation is only valid for incompressible flows (ie.
only if div(U)=0). It is obtained by expanding the conservative formulation
and simplifying with div(U)=0. The linear combination of gradients comes from
this first equation and the fact that (grad(U)-grad(U)^T).W = 0.
"""
check_instance(formulation, StretchingFormulation)
check_instance(velocity, Field)
check_instance(vorticity, Field)
check_instance(dt, ScalarParameter)
check_instance(variables, dict, keys=Field, values=CartesianTopologyDescriptors)
check_instance(base_kwds, dict, keys=str, allow_none=True)
check_instance(name, str, allow_none=True)
if (velocity.dim != 3) or (velocity.nb_components != 3):
msg = "Given velocity is not a 3d vector field."
raise ValueError(msg)
if (vorticity.dim != 3) or (vorticity.nb_components != 3):
msg = "Given vorticity is not a 3d vector field."
raise ValueError(msg)
name = first_not_None(name, "stretching")
base_kwds = first_not_None(base_kwds, {})
C = first_not_None(C, sm.Integer(1))
if isinstance(C, npw.ndarray):
assert C.ndim in (0, 1), C.ndim
assert C.size in (1, 3), C.size
if formulation is StretchingFormulation.MIXED_GRAD_UW:
A = first_not_None(A, sm.Rational(1, 2))
if isinstance(A, npw.ndarray):
assert A.ndim == 0, A.ndim
assert A.size == 1, A.size
elif A is not None:
msg = f"Formulation is {formulation} but A was specified."
raise ValueError(msg)
elif formulation is StretchingFormulation.GRAD_UW:
A = sm.Integer(1)
formulation = StretchingFormulation.MIXED_GRAD_UW
elif formulation is StretchingFormulation.GRAD_UW_T:
A = sm.Integer(0)
formulation = StretchingFormulation.MIXED_GRAD_UW
exprs = self._gen_expressions(formulation, velocity, vorticity, C, A)
super().__init__(
name=name,
variables=variables,
dt=dt,
candidate_input_tensors=(
velocity,
vorticity,
),
candidate_output_tensors=(vorticity,),
base_kwds=base_kwds,
implementation=implementation,
exprs=exprs,
**kwds,
)
from hysop.operator.base.custom_symbolic_operator import (
CustomSymbolicOperatorBase,
)
if not issubclass(self._operator, CustomSymbolicOperatorBase):
msg = "Class {} does not inherit from the directional operator interface "
msg += "({})."
msg = msg.format(self._operator, CustomSymbolicOperatorBase)
raise TypeError(msg)
def _gen_expressions(self, formulation, velocity, vorticity, C, A):
from hysop.symbolic.base import TensorBase
from hysop.symbolic.field import div, grad
frame = velocity.domain.frame
U = velocity.s()
W = vorticity.s()
lhs = W.diff(frame.time)
if formulation is StretchingFormulation.CONSERVATIVE:
assert A is None, A
UW = npw.outer(U, W).view(TensorBase).freeze()
# Freezing is important because sympy will split the derivatives,
# ie. makes the stretching term not conservative anymore
# (and time integration *WILL* fail)
# See Taylor-Green Re=1600 without calling freeze().
rhs = div(UW, frame)
elif formulation is StretchingFormulation.MIXED_GRAD_UW:
assert A is not None, A
rhs = (A * grad(U, frame) + (sm.Integer(1) - A) * grad(U, frame).T).dot(W)
else:
msg = f"Unknown stretching formulation {formulation}."
raise ValueError(msg)
rhs = C * rhs
exprs = Assignment.assign(lhs, rhs)
return exprs
[docs]
class StaticDirectionalStretching(DirectionalOperatorFrontend):
"""
Directional stretching using the symbolic code generation framework.
"""
[docs]
@classmethod
def implementations(cls):
from hysop.backend.host.python.operator.directional.stretching_dir import (
PythonDirectionalStretching,
)
__implementations = {Implementation.PYTHON: PythonDirectionalStretching}
return __implementations
[docs]
@classmethod
def default_implementation(cls):
return Implementation.PYTHON
@debug
def __new__(
cls,
formulation,
velocity,
vorticity,
variables,
dt,
C=None,
A=None,
name=None,
implementation=None,
base_kwds=None,
**kwds,
):
base_kwds = first_not_None(base_kwds, {})
return super().__new__(
cls,
name=name,
variables=variables,
dt=dt,
formulation=formulation,
velocity=velocity,
vorticity=vorticity,
C=C,
A=A,
candidate_input_tensors=None,
candidate_output_tensors=None,
base_kwds=base_kwds,
implementation=implementation,
**kwds,
)
@debug
def __init__(
self,
formulation,
velocity,
vorticity,
variables,
dt,
C=None,
A=None,
name=None,
implementation=None,
base_kwds=None,
**kwds,
):
"""
Initialize directional stretching operator frontend.
Vortex stretching is the lengthening of vortices in three-dimensional fluid flow,
associated with a corresponding increase of the component of vorticity in the stretching
direction due to the conservation of angular momentum.
Solves dW/dt = C * f(U,W)
where U is the velocity, W the vorticity and f is one of the following formulation:
CONSERVATIVE FORMULATION: f(U,W) = div( U : W )
MIXED GRADIENT FORMULATION: f(U,W) = (A*grad(U) + (1-A)*grad(U)^T) . W
where : represents the outer product U:W = (Ui*Wj)ij.
. represents the matrix-vector product.
^T is the transposition operator
U and W are always three dimensional fields.
C is scalar or 3d vector-like of symbolic coefficients.
A is a scalar symbolic coefficient.
f(U,W) is always directionally splittable.
Parameters
----------
formulation: hysop.constants.StretchingFormulation
The formulation of this stretching operator:
CONSERVATIVE: f(U,W) = div( U : W )
GRAD_UW: f(U,W) = grad(U) . W
GRAD_UW_T: f(U,W) = grad(U)^T . W
MIXED_GRAD_UW: f(U,W) = (A*grad(U) + (1-A)*grad(U)^T) . W
where A is a user given scalar or 3d vector like that contains
scalar (floating point) coefficient(s) or parameter(s).
The most common mixed gradient formulation is A=0.5.
GRAD_UW, GRAD_UW_T and MIXED_GRAD_UW formulations are only valid
for incompressible flows, see Notes.
velocity: Field
Velocity U as read-only input three-dimensional continuous field.
vorticity: Field
Vorticity W as read-write three-dimensional continuous field.
variables: dict
Dictionary of fields as keys and topology descriptors as values.
dt: ScalarParameter
Timestep parameter that will be used for time integration.
C: scalar or vector like of 3 symbolic coefficients, optional
The stretching leading coefficient C can be scalar or vector like.
Contained values should be numerical coefficients, parameters or
generic symbolic expressions such that C*f(U,W) is directionally splittable.
Here * is the classical elementwise multiplication.
Default value is 1.
A: scalar symbolic coefficient, optional
Should only be given for MIXED_GRAD_UW formulations.
ValueError will be raised on other formulations.
The linear combination coefficients A is a scalar.
Contained value should be a numerical coefficient, a parameter (or a generic symbolic
expression) such that C*(A*grad(U).W + (1-A)*grad^T(U).W) is directionally splittable.
Here * is the classical elementwise multiplication.
Default value is 0.5.
name: str, optional, defaults to 'stretching'.
Name of this stretching operator.
implementation: Implementation, optional, defaults to None
Target implementation, should be contained in available_implementations().
If None, implementation will be set to default_implementation().
base_kwds: dict, optional, defaults to None
Base class keywords arguments.
If None, an empty dict will be passed.
kwds:
Keywords arguments that will be passed towards implementation
operator __init__.
Notes
-----
The MIXED GRADIENT formulation is only valid for incompressible flows (ie.
only if div(U)=0). It is obtained by expanding the conservative formulation
and simplifying with div(U)=0. The linear combination of gradients comes from
this first equation and the fact that (grad(U)-grad(U)^T).W = 0.
"""
check_instance(formulation, StretchingFormulation)
check_instance(velocity, Field)
check_instance(vorticity, Field)
check_instance(dt, ScalarParameter)
check_instance(variables, dict, keys=Field, values=CartesianTopologyDescriptors)
check_instance(base_kwds, dict, keys=str, allow_none=True)
check_instance(name, str, allow_none=True)
if (velocity.dim != 3) or (velocity.nb_components != 3):
msg = "Given velocity is not a 3d vector field."
raise ValueError(msg)
if (vorticity.dim != 3) or (vorticity.nb_components != 3):
msg = "Given vorticity is not a 3d vector field."
raise ValueError(msg)
name = first_not_None(name, "stretching")
base_kwds = first_not_None(base_kwds, {})
C = first_not_None(C, sm.Integer(1))
if isinstance(C, npw.ndarray):
assert C.ndim in (0, 1), C.ndim
assert C.size in (1, 3), C.size
if formulation is StretchingFormulation.MIXED_GRAD_UW:
A = first_not_None(A, sm.Rational(1, 2))
if isinstance(A, npw.ndarray):
assert A.ndim == 0, A.ndim
assert A.size == 1, A.size
elif A is not None:
msg = f"Formulation is {formulation} but A was specified."
raise ValueError(msg)
elif formulation is StretchingFormulation.GRAD_UW:
A = sm.Integer(1)
formulation = StretchingFormulation.MIXED_GRAD_UW
elif formulation is StretchingFormulation.GRAD_UW_T:
A = sm.Integer(0)
formulation = StretchingFormulation.MIXED_GRAD_UW
super().__init__(
name=name,
variables=variables,
dt=dt,
formulation=formulation,
velocity=velocity,
vorticity=vorticity,
C=C,
A=A,
candidate_input_tensors=(
velocity,
vorticity,
),
candidate_output_tensors=(vorticity,),
base_kwds=base_kwds,
implementation=implementation,
**kwds,
)